Two-step melting in two dimensions: First-order liquid- hexatic transition 
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Melting in two spatial dimensions, as realized in thin films or at interfaces, represents one of 
the most fascinating phase transitions in nature, but it remains poorly understood. Even for the 
fundamental hard-disk model, the melting mechanism has not been agreed on after fifty years of 
studies. A recent Monte Carlo algorithm allows us to thermalize systems large enough to access the 
thermodynamic regime. We show that melting in hard disks proceeds in two steps with a liquid 
phase, a hexatic phase, and a solid. The hexatic-solid transition is continuous while, surprisingly, 
the liquid-hexatic transition is of first-order. This melting scenario solves one of the fundamental 
statistical-physics models, which is at the root of a large body of theoretical, computational and 
experimental research. 



Generic two-dimensional particle systems cannot crys- 
tallize at finite temperature [lj— [3(] because of the impor- 
tance of fluctuations, yet they may form solids This 
paradox has provided the motivation for elucidating the 
fundamental melting transition in two spatial dimen- 
sions. A crystal is characterized by particle positions 
which fluctuate about the sites of an infinite regular lat- 
tice. It has long-range positional order. Bond orienta- 
tions are also the same throughout the lattice. A crystal 
thus possesses long-range orientational order. The po- 
sitional correlations of a two-dimensional solid decay to 
zero as a power law at large distances. Because of the 
absence of a scale, one speaks of "quasi-long range" or- 
der. In a two-dimensional solid, the lattice distortions 
preserve long-range orientational order [5], while in a liq- 
uid, both the positional and the orientational correlations 
decay exponentially. 

Besides the solid and the liquid, a third phase, called 
"hexatic" , has been discussed but never clearly identified 
in particle systems. The hexatic phase is characterized by 
exponential positional but quasi-long range orientational 
correlations. It has long been discussed whether the melt- 
ing transition follows a one-step first-order scenario be- 
tween the liquid and the solid (without the hexatic) as in 
three spatial dimensions [6]), or whether it agrees with the 
celebrated Kosterlitz, ThoulessQ, Halperin, Nelsonjij 
and Young [9] (KTHNY) two-step scenario with a hex- 
atic phase separated by continuous transitions from the 
liquid and the solid [13413. 

Two-dimensional melting was discovered Q in the sim- 
plest particle system, the hard-disk model. Hard disks (of 
radius a) are structureless and all configurations of non- 
overlapping disks have zero potential energy. Two iso- 
lated disks only feel the hard-core repulsion, but the other 
disks mediate an entropic "depletion" interaction (see, 
e.g., [ll|). Phase transitions result from an "order from 
disorder" phenomenon: At high density, ordered configu- 
rations can allow for larger local fluctuations, thus higher 
entropy, than the disordered liquid. For hard disks, no 



difference exists between the liquid and the gas. At fixed 
density 77, the phase diagram is independent of tempera- 
ture T = 1/&bA and the pressure is proportional to T, 
as discovered by D. Bernoulli in 1738. Even for this basic 
model, the nature of the melting transition has not been 
agreed on. 

The hard-disk model has been simulated with the 
local Monte Carlo algorithm since the original work 
by Metropolis et al. [20|. A faster collective- move 
"event-chain" Monte Carlo algorithm was developed only 
recently (see [22|). We will use it to show that the 
melting transition neither follows the one-step first-order 
nor the two-step continuous KTHNY scenario. 

To quantify orientational order, we express the local 
orientation of disk k through the complex vector = 
(exp(6z0fcz)), with () the average over all the neighbors I 
of k. The angle fai describes the orientation of the bond 
kl with respect to a fixed axis. The sample orientation 
is defined as \£ = l/N^2 k S&k- For a perfect triangular 
lattice, all the angles 6(j)ki are the same and = = 
1 (see [22[). 

In Fig. [TJ the local orientations of a configuration with 
N = 1024 2 disks at density n = Nna 2 /V = 0.708 in 
a square box of volume V are projected onto the sam- 
ple orientation and represented using a color code (see 
[22[). Inside this configuration, a vertical stripe with 
density ~ 0.716 preserving orientational order over long 
distances coexists with a stripe of disordered liquid of 
lower density ~ 0.700. Each stripe corresponds to a dif- 
ferent phase. The two interfaces of length c± \/N close on 
themselves via the periodic boundary conditions. Stripe- 
shaped phases as in Fig. [T^i are found in the center of a 
coexistence interval n G [0.700,0.716], whereas close to 
its endpoints, a "bubble" of the minority phase is present 
inside the majority phase for n > 0.700 and n < 0.716 
(see Fig. [2]). This phase coexistence is the hallmark of a 
first-order transition. 

The first-order transition shows up in the equilib- 
rium equation of state P(V) (see Fig. [2]). At finite N, 
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FIG. 1. Phase-coexistence for 1024 2 thermalized hard disks at density 77 = 0.708. a: Color-coded local orientations 
showing long orient at ional correlations (blue region, see b,c) coexisting with short-range correlations (see dV e: Local densities 
(averaged over a radius of 50cr), demonstrating the connection between density and local orientation (see [22j]). In b, c, d, disks 
with five (seven) neighbors are colored in gray (black). 



the free energy is not necessarily convex (as it would 
be in an infinite system) and the equilibrium pressure 
P(V) = —dF/dV can form a thermodynamically stable 
loop due to the interface free energy. The pressure loop 
in the coexistence window of a finite system is caused by 
the curved interface between a bubble of minority phase 
and the surrounding majority phase (see Fig. |2)b,d)). In 
a system with periodic boundary conditions, the pres- 
sure loop contains a horizontal piece corresponding to 
the "stripe" regime, where the interfaces are flat. This is 
visible near 77 ~ 0.708 for the largest systems in Fig. [2j In 
a finite system, the Maxwell construction suppresses the 
interface effects. For the equation of state of Fig. [2^, this 
construction confirms the boundary densities 77 = 0.700 
and rj — 0.716 of Fig. [T] for the coexistence interval, with 
very small finite-size effects. The interface free energy per 
disk, the hatched area in Fig. [2j depends on the length 
oc VN of the interface in the "stripe" regime so that 
A/ = AF/N oc 1/y/N (see Fig. EF). 

The first-order nature of the transition involving the 
liquid is thus established by i): The visual evidence of 
phase coexistence in Fig. [TJ ii): The oc 1/y/N scaling 
of the interface free energy per disk(23|. and Hi): The 
characteristic shape of the equation of state in a finite 
periodic system [2j426|. We stress that the system size 
is larger than the physical length scales so that the results 
hold in the thermodynamic limit (see |22j). 

In the coexistence interval, the individual phases are 
difficult to analyze at large length scales because of the 
fluctuating interface, and only the low-density coexisting 



phase is identified as a liquid with orientational corre- 
lations below a scale of ~ 100a (see Fig. QJi,d). Un- 
like constant NV simulations, Gibbs ensemble simula- 
tions can have phase coexistence without interfaces, but 
these simulations are very slow at large TV (see [22j]). The 
single-phase system at density 77 = 0.718, is above the co- 
existence window for all N (see Fig. [2]), and it allows us 
to characterize the high-density coexisting phase. 

Positional order can be studied in the two-dimensional 
pair correlation g(Ar), the high-resolution histogram of 
periodic pair distances Ar^ = — Vj sampled from all 
N(N — l)/2 pairs i,j of disks. To average this two- 
dimensional histogram over configurations (as in Fig. [3]) 
the latter are oriented such that the Ax axis points in the 
direction of the sample orientation At short distances, 
hexagonal order is evident at 77 = 0.718 (see Fig.[3k). The 
excellent contrast between peaks and valleys of g(Ar) at 
small I Ar| > 2a underlines the single-phase nature of the 
system at this density. The cut of the histogram along the 
positive Ax axis leaves no doubt that the system has ex- 
ponentially decaying positional order on a length scale of 
^ lOOcr and cannot be a solid. The (one-dimensional) po- 
sitional correlation function Ck(r), computed by Fourier 
transform of g(Ar), fully confirms these statements (see 

0). 

The orientational correlations at density 77 = 0.718 de- 
cay extremely slowly and do not allow us to distinguish 
between quasi- long range and long-range order (see [22|). 
However, short-ranged positional correlation is inconsis- 
tent with long-ranged orientational order. It follows that 
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FIG. 2. Equilibrium equation of state for hard disks. The 
pressure is plotted vs. volume per particle (v = V/N) (lower 
scale) and density 77 (upper scale)). In the coexistence region, 
the strong system-size dependence stems from the interface 
free energy. The Maxwell constructions (horizontal lines) sup- 
press the interface effects (with a convex free energy) for each 
N. "Stripe" (c, for N = 1024 2 ) and "bubble" configurations 
(b,d) are shown in the coexistence region, together with two 
single-phase configurations (a,e). The interface free energy 
per disk f3Af (hatched area) scales as 1/y/N (f). 



the orientation must be quasi-long ranged with a small 
exponent < 0, and that the system at 77 = 0.718 and the 
high-density coexisting phase are both hexatic. 

The two-dimensional pair correlation g(Ar) — 1 of 
Fig. 03 allows us to follow the transition from the hexatic 
to the solid: The positional order increases continuously 
with density and crosses over into power-law behavior at 
density 77 ~ 0.720, with an exponent ~ — 1/3 which cor- 
responds to the stability limit of the solid phase in the 
KTHNY scenario. The hexatic-solid transition thus takes 
place at 77 > 0.720. At this density, the positional correla- 
tion function at large distances r, displays the finite-size 
effects characteristic of a continuous transition, but up 
to a few hundred a, Ck is well stabilized with system size 
(see [Hj]). Moreover, no pressure loop is observed in the 
equation of state, and the compressibility remains very 
small. The system is clearly in a single phase. Unlike 
the liquid-hexatic transition, the hexatic-solid transition 
therefore follows the KTHNY scenario, and is continu- 



ous. 

The single-phase hexatic regime is confined to a den- 
sity interval 77 G [0.716,0.720]. Although narrow, it is an 
order of magnitude larger than the scale set by density 
fluctuations for our largest systems and can be be easily 
resolved (see (Hj). In the hexatic phase, the orientational 
correlations decay extremely slowly. The exponent of the 
orientational correlations is close to zero and negative. It 
remains far from the lower limit of —1/4 at the contin- 
uous KTHNY transition, as this transition is preempted 
by a first-order instability. 

The event-chain algorithm is about two orders of mag- 
nitude faster than the local Monte Carlo used up to now, 
allowing us to thermalize for the first time dense sys- 
tems with up to 1024 2 disks. To illustrate convergence 
toward thermal equilibrium and to check that hard disks 
in the window of densities 77 G [0.700,0.716] are indeed 
phase-separated, we show in Fig. |4] two one- week simula- 
tions of our largest systems after quenches from radically 
different initial conditions, namely the (unstable) crys- 
tal, with |^| = 1, and the liquid, for which |^| ~ 0. For 
both initial conditions, a slow process of coarsening takes 
place (see Fig. Hk,b). Phase separation is observed after 
~ 10 6 displacements per disk, and the sample orienta- 
tion takes on similar absolute values (see Fig. |4fc). Ef- 
fective simulation times of many earlier calculations were 
much shorter 14 , li| , and the simulations remained in an 
out-of-equilibrium state which is homogeneous on large 
length scales, whereas the thermalized system is phase- 
separated and therefore inhomogeneous. The production 
runs for N = 1024 2 were obtained from Markov chains 
with running times of nine months, 30 times larger than 
those of Fig. Hk,b. 

The solution of the melting problem presented in this 
work provides the starting point for the understanding of 
melting in films, suspensions, and other soft-condensed- 
matter systems. The insights obtained combine thermo- 
dynamic reasoning with powerful tools: advanced simula- 
tion algorithms, direct visualization, and a failsafe anal- 
ysis of correlations. These tools will all be widely ap- 
plicable, for example to study the cross-over from two to 
three-dimensional melting as it is realized experimentally 
with spheres under different confinement conditions [I?}- 

In simple systems such as hard disks and spheres, en- 
tropic and elastic effects have the same origin: elastic 
forces are entropically induced. For general interaction 
potentials, entropy and elasticity are no longer strictly 
linked and order-disorder transitions, which can then 
take place as a function of temperature or of density, 
might realize other melting scenarios [27]. Theoretical, 
computational and experimental research on more com- 
plex microscopic models will build on the hard-disk so- 
lution obtained in this work. 

We are indebted to K. Binder and D. R. Nelson for 
helpful discussions and correspondence. We thank J. 
Dalibard and G. Bastard for a critical reading of the 




FIG. 3. Configuration-averaged two-dimensional pair correlation. g(Ar) is obtained from the two-dimensional histogram 
of periodic distances Ar ij = n — Yj. a: Pair correlation g(Ar) at density 77 = 0.718 for small Ar = (Ax, Ay). Each disk 
configuration is oriented with respect to \P. The excellent contrast between the peak and the bottom values of g(Ar) at 
|Ar| > 2cr, of about (16 : 0.2), provides evidence for the single-phase nature of the system, b: Cut of the sample- averaged 
g(Ar) — 1 for Ar = (Ax, 0). Decay is exponential for 77 = 0.718 and algebraic for 77 = 0.720. (See [52] for positional and 
orient at ional correlation functions.) 




FIG. 4. Approach to thermal equilibrium from different initial conditions. a,b: 1024 2 hard disks at density 77 = 0.708, after 
a quench from a high-density crystal (a) and from a low-density liquid (b), showing coarsening leading to phase separation 
(Color code for as in Fig. [TJd, see also [22j]). Each of the runs takes about one week of CPU time, c: Absolute value of the 
sample orientation for the simulations in a,b, compared to runs with the local Monte Carlo algorithm from the same initial 
conditions (time in attempted displacements per disk). The correlation time of the event-chain algorithm, on the order of 10 6 
displacements per disk, estimated from c, agrees with the correlation time estimated in our production runs with 6 x 10 7 total 
displacements per disk. 
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